REPORT  DOCUMENTATION  PAGE 

Form  Approved  0MB  No.  0704-0188 

Public  reporting  burden  for  this  collection  of  infomnaiion  is  estimated  to  average  1  hour  per  resp 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  inforr 
collection  of  information,  including  suggestions  for  reducing  this  burden  to  Washington  Headqua 
navis  Hiohwav.  Suite  1204.  Arlinaton.  VA  22202-4302,  and  to  the  Office  of  Management  and  Bud( 

snse,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
nation.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
rters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson 
jet.  Paperwork  Reduction  Project  (0704-0188),  Washington,  DC  20503. 

1 .  AGENCY  USE  ONLY  (Leave  blank) 

2.  REPORT  DATE 

3  December  1996 

3.  REPORT  TYPE  AND  DATES  COVERED 

Final  Report 

4.  TITLE  AND  SUBTITLE 


Computational  Damage  Model  for  Layered  Composite  Materials  (Numerical  Modeling  of  Dynamical 
Deforming  of  Damageable  Thermoviscoelastic  Composite  Shell  in  Internal  Loading) 


6.  AUTHOR(S) 

Prof.  Nickoiay  Smirnov 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Moscow  State  University 
Subdept  of  Gas  and  Wave  Dynamics 
Moscow  1 1 9899 
Russia 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 
EOARD 

PSC  802  BOX  14 
FPO  09499-0200 


F6170896W0103 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

N/A 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 

SPC  96-4012 


11.  SUPPLEMENTARY  NOTES 


12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 


12b.  DISTRIBUTION  CODE 
A 


13.  ABSTRACT  (Maximum  200  words) 

This  report  results  from  a  contract  tasking  Moscow  State  University  as  follows:  The  contractor  will  create  a  new  mathematical  model  for 
dynamical  deforming  and  failure  of  viscoelastic  damageable  unidirectional  laminated  composites  as  detailed  in  his  proposal  dated  28  Dec 
95. 


19910113  005 


14.  SUBJECT  TERMS 


Aircraft  Subsystem 


15.  NUMBER  OF  PAGES 
165 


16.  PRICE  CODE 
N/A 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

UNCLASSIFIED 


1 8.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

UNCLASSIFIED 


19,  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

UNCLASSIFIED 


20.  LIMITATION  OF  ABSTRACT 


UL 


NSN  7540-01-280-5500 


DTic  q;: 


.llZD  1 


Standard  Form  298  (Rev.  2-89) 
Prescribed  by  ANSI  Std.  239-18 
298-102 


r 


TfflS  DOCUMENT  IS  BEST 
QUALITY  AVAILABLE.  THE  COPY 
FURNISHED  TO  DTIC  CONTAINED 
A  SIGNIFICANT  NUMBER  OF 
COLOR  PAGES  WHICH  DO  NOT 
REPRODUCE  LEGIBLY  ON  BLACK 
AND  WHITE  MICROnCHE. 


Final  Report 
on 

SPC  96-4012 

COMPUTATIONAL  DAMAGE  MODEL  FOR 
LAMINATED  COMPOSITE  MATERIALS 

Numerical  Modelling  of  Dynamical  Deforming  of  Damageable 
Thermoviscoelastic  Composite  Shell  in 
Internal  Loading. 


Principal  Investigator: 
Prof.,Dr.Sc.(hab.)  N.N. Smirnov 


Investigators: 

M.Sc.  V.F. Nikitin 
Prof.,Dr.Sc.(hab.)  A.B.Kiselev 
Dr.  I .D  .Dimitrienko 
Dr.  V.R.Dushin 


19970113  005 

BRUSSELS  -  1996 


DTIC  QUALmr  mSPECl'SD  i 


Contents 


Introduction .  2 

1  Part  I.  Mathematical  model  for  damageable  thermoviscoelastic  lamina¬ 
ted  composite  shell.  4 

1.1  Foreword .  4 

1.2  Governing  system  of  equations  for  dynamical  deforming  of  an 

axisymmetrical  composite  containment .  5 

1.3  Numerical  modelling .  11 

2  Part  II.  Results  of  numerical  modelling.  14 

2.1  Foreword .  14 

2.2  Determining  the  internal  loading  of  a  containment .  15 

2.2.1  Mass  balance  equations  in  the  gas  phase .  15 

2.2.2  Momentum  equation  for  the  gas  phase .  16 

2.2.3  Energy  balance  equation  for  the  gas  phase .  17 

2.2.4  Fluxes  modelling  and  closing  the  K-epsilon  model .  18 

2.2.5  Initial  and  boundary  conditions .  19 

2.2.6  Computational  techniques .  19 

2.2.7  Results  of  numerical  modelling  of  internal  loading .  21 

2.3  Dynamical  deforming  of  walls  and  accumulation  of  damages.  .  .  54 

2.3.1  Evolution  of  the  shell  parameters  in  uniform  dynamical 

loading .  54 

2.3.2  Shell  parameters  in  nonuniform  dynamical  loading .  113 

References .  163 

Conclusions .  164 


1 


INTRODUCTION 


The  results  of  theoretical  investigations  undertaken  within  the  frames  of  the  present 
project  brought  to  life  a  new  approach  to  modelling  of  thermoviscoelastic  damageable 
composite  materials  with  laminated  or  fibrous  structures.  The  mathematical  models  for 
dynamical  deforming  and  failure  of  such  composite  materials  is  described  in  details  in  the 
interim  reports  [1,2],  The  main  results  of  the  investigations  are  the  following. 

The  structural  breakup  of  composite  materials  has  three  characteristic  stages: 

1.  Loading  of  the  structure  and  accumulation  of  irreversible  damages.  2.  Fracture 
occuring  in  the  zones  where  failure  criterion  is  satisfied.  3.  Post-fracture  deformations, 
formation  of  cracks  and  fragments. 

At  the  present  phase  of  the  Project  the  main  attention  was  paid  to  the  first  two 
stages  of  the  process.  It  was  stated  that  failure  generated  from  quasi-statically  applied 
over-pressures  differs  significantly  from  dynamically  applied  over-pressures.  Quasi-static 
matrix  cracking  may  lead  to  ’’failure”  in  the  form  of  a  leak  or,  as  a  worst  case,  catastrophic 
failure  involving  matrix  cracking  and  formation  of  segments  attached  to  each  other  by 
the  fibers.  Pressures  applied  dynamically  (with  relatively  high  amplitudes  and  rates  of 
loadings)  force  multiple  fractures  in  all  of  the  wall’s  constituent  materials  and  result  in 
the  generation  of  many  small  particles. 

A  two-phase  failure  model  was  worked  out  for  composite  materials  wherein  irreversible 
work  for  damage  production  was  used  to  determine  the  extent  of  damage  accumulated 
in  local  regions  of  the  structure.  The  model  computes  irreversible  damage  based  on  the 
growth  of  pores  and  the  formation  of  dislocations  and  delaminations  of  layers.  The  cha¬ 
racteristics  of  these  processes  are  modeled  as  scalar  quantities.  The  amount  of  damage 
that  a  material  experiences  is  measured  by  the  specific  dissipation  D.  All  of  the  indi¬ 
vidual  damage  processes  are  lumped  together  into  the  global  quantity  D.  Failure  of  a 
material  is  determined  by  comparing  D  to  an  experimentally  determined  material  pro¬ 
perty  D* .  D*  represents  the  summation  of  all  irrecoverable  material  motions  upon  failure 
of  the  laminate.  D*  is  presently  a  volummetric  equivalent  failure  criterion  analogous  to  a 
spherical  failure  envelope.  The  value  of  D*  remaines  unknown  in  the  equations  [1,2]  but 
experimental  guidance  describing  one  of  the  possible  methods  of  determining  the  scalar 
value  of  D*  for  spherical  envelope  was  provided  in  the  Third  Quaterly  Report  [3]. 

The  worked  out  failure  model  considers  3-D  stresses.  Specific  dissipation  D  is  being 
accumulated  in  each  space  element  as  the  pressure  loading  continues.  When  D  =  D*  is 
attained  the  element  is  assumed  to  have  locally  failed.  The  zones  of  the  structure  wherein 
the  failure  criterion  D  —  D*  is  satisfied  are  considered  to  be  damaged.  A  detailed  fracture 
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progression  model  and  methods  of  determining  of  fragments’  characteristics  have  not  been 
worked  out  yet  for  the  thermoviscoelastic  laminated  composite  materials. 

To  demonstrate  the  utility  of  the  model  a  problem  of  dynamical  deforming  and  for¬ 
mation  of  the  failure  zones  in  an  axisymmetrical  shell  of  the  laminate  under  the  influence 
of  internal  loading  was  regarded. 

The  governing  system  of  equations  for  thermoviscoelastic  laminates  was  written  in  the 
lagrangian  variables  to  be  applied  for  the  description  of  behaviour  of  an  axisymmetrical 
shell.  The  numerical  model  describing  the  dynamical  deforming  of  the  shell  made  use  of 
a  finite  element  approach.  To  determine  the  rates  of  internal  wall  loading  the  model  was 
incorporated  into  a  hydrocode  enabling  to  solve  the  problems  of  shock  waves  propagation 
and  reflections  in  gas  in  a  closed  vessel  after  internal  explosion. 

The  Part  I  describes  the  mathematical  model  being  used  for  the  solution  of  the  problem 
of  dynamical  deforming  of  an  axisymmetrical  laminated  composite  containment  under  the 
influence  of  internal  loading  in  reflected  shock  waves.  The  model  is  based  on  the  assump¬ 
tions  of  the  thin  shell  theory  and  enables  to  evaluate  the  accumulation  of  damges  in  the 
shell  and  the  growth  of  spesific  dissipation  and  damage  parameters  responsible  for  dama¬ 
ges  in  tension,  shear  and  delamination.  The  model  incorporates  three  additional  material 
constants  characterizing  the  resistance  of  the  laminate  to  accumulating  the  damages.  The 
numerical  method  is  also  described  briethly. 

The  Part  II  describes  the  results  of  numerical  modelling  of  several  problems  characte¬ 
rised  by  different  internal  loading  pattern.  The  numerical  solution  of  those  problems  was 
undertaken  to  demonstrate  the  utility  of  the  worked  out  model  for  damageable  laminated 
composite  material  being  incorporated  into  a  hydrocode. 
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Chapter  1 

Part  I.  Mathematical  model  for 
damageable  thermoviscoelastic 
laminated  composite  shell. 


1.1  Foreword. 


The  pjresent  Part  describes  the  model  for  damageable  thermoviscoelastic  laminated 
composite  material  adapted  for  solving  the  problems  of  dynamical  deforming  and  failure 
of  axisymmetrical  thin-walled  containments  in  internal  loading.  Two  possible  cases  of 
dynamical  loading  are  regarded:  uniform  and  nonuniform  loading.  In  the  last  case  the 
loads  are  determined  with  the  help  of  a  hydrocode  modelling  shock  waves  propagation  in 
gas  inside  a  closed  vessel  and  reflections  from  the  walls  of  the  vessel.  The  shock  waves 
inside  the  vessel  originate  due  to  a  definite  discontinuity  in  the  initial  conditions  simulating 
an  explosion.  The  initial  zone  of  high  pressure  and  density  gas  can  be  located  in  any  place 
on  the  axis  of  symmetry  and  be  of  a  variable  volume. 

The  dynamical  loads  bring  to  an  expansion  of  the  walls  of  the  containment  and  accu¬ 
mulations  of  damages.  Since  the  actual  criterion  (the  value  of  critical  specific  dissipation 
D*  for  composites  remains  unknown  the  numerical  programme  continues  working  after 
the  assumed  critical  values  of  dissipation  are  surpassed  in  some  zones  or  even  everywhere. 
To  terminate  the  programme  until  the  results  loose  the  physical  sence  additional  criteria 
are  used  making  it  possible  to  avoid  overshoot  in  the  damage  parameters.  That  makes  it 
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possible  for  us  to  follow  the  dynamics  of  the  evolution  of  the  zones  in  the  shell  wherein  the 
values  of  the  specific  dissipation  surpass  the  assumed  critical  values.  In  fact  the  nume¬ 
rical  simulations  based  on  the  present  model  should  be  terminated  as  soon  as  D  reaches 
critical  value  for  anyone  of  the  elements  because  the  model  does  not  take  into  account 
the  changes  in  the  stresses  and  strains  of  the  neighbouring  elements  in  case  of  a  failure  of 
one  of  the  elements.  The  results  of  the  present  modelling  are  destined  to  demonstrate  the 
applicability  of  the  model  and  to  serve  the  base  for  the  further  development  of  composite 
breakup  modelling. 


1.2  Governing  system  of  equations  for  dynamical 
deforming  of  an  axisymmetrical  composite  con¬ 
tainment. 


The  system  of  equation  derived  in  the  present  chapter  was  obtained  making  use  of  the 
model  and  notations  described  in  details  in  [2].  Thus  the  general  description  of  the  model 
is  not  given  here.  If  one  applies  results  of  our  previous  investigations  (see  reports  [1-3]) 
to  an  axisymmetrical  thick  shell,  one  will  obtain  the  following  results. 

Let  x^r  be  the  physical  cylindrical  coordinate  system,  and  z,s  -  the  lagrangian  coor¬ 
dinates  attached  to  the  shell,  so  that  the  coordinate  s  changes  along  the  generating  line 
of  the  shell  and  2  -  across  the  shell.  The  equation  5;  =  0  is  satisfied  on  the  surface  gene¬ 
rated  by  a  set  of  points  laying  in  the  middle  of  the  shell  thickness.  The  part  2  >  0  faces 
the  external  side  of  the  shell  (see  Fig.l.)  The  middle  surface  of  the  shell  is  represented 
by  function  r  =  ro{x)  in  the  physical  coordinate  system,  or  the  parametrical  system  of 
equations; 

r  =  ro(5),  X  =-  Xo{s). 

It  is  assumed  that  the  displacement  of  the  shell  is  represented  in  the  following  terms: 


u{z,  s,  t)  =  w°(s,  t)  —  zG(s,  t),  w(z,  s,  t)  =  w°{s,  t), 


where  u  is  tangential  displacement  directed  along  the  generating  line,  and  w  is  the  normal 
displacement  directed  collinear  to  the  external  normale  to  the  shell. 

With  this  assumption,  the  system  of  motion  equations  gets  the  following  form: 


_  1  djroNi) 

^  ?’o  ds 


-^N2  +  QKi 
ro  ds 
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-p 


d'^vP 

Vl'W 


1  ^(?oQ) 
ro  9s 
1  djroMi 
Tq  ds 


-  NiK,  -  N2K2  +  -^Pni-h/2) 

- —M2  -  Q 

ro  os 


Here  p  =<  p  >—  cpi  +  (1  —  c)p2  is  the  effective  density  of  material  and  the  other  terms 
in  the  right  side  are  the  following: 

Ki,K2  are  the  physical  main  curvatures  of  the  shell’s  middle  surface,  so  that  Ki  is 
inverse  to  the  radius  of  curvature  along  the  generating  line,  and  K2  is  the  curvature 
in  the  perpendicular  direction.  The  curvatures  are  determined  in  terms  of  the  physical 
coordinates  ?’°(s,t),  x^{s,t)  of  the  middle  surface: 


Ki 


drg  d'^XQ  dxo  d-rg 

ds  ds~  ds  ds'^ 


Ni,  N2  are  the  shell  tensions  in  the  directions  along  the  generating  line  and  along  the 
orthogonal  geodesical  line  correspondingly: 


h/2 


an  dz, 


-h/2 


where  h  is  the  initial  thickness  of  the  shell; 

Q  is  the  rotational  tension  determined  by; 


-h/2 


Ml,  M2  are  the  bending  moments  along  the  generating  line  and  perpendicular  geode¬ 
sical  line  correspondingly: 

1 


M  = 


h/2 

M  f 

h  J 


zan  dz. 


-h/2 

These  tensions  and  moments  are  determined  in  terms  of  elastic  deformations  tensor 
components,  enthropy  and  damadge  parameters.  But  the  elastic  deformations  themselves 
are  determined  by  the  values  of  full  deformations  tensor  components. 
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With  our  assumption  of  linear  distribution  of  functions  across  the  shell  thickness,  the 
components  of  the  full  deformations  tensor  velocities  are  represented  in  the  following  way; 

ei{z,s,t)  =  -  zki(s,t),  €2(^,5, t)  =  “  zk2{s,t),  iu  = 

The  linkage  between  the  velocities  of  displacements  and  full  deformations’  velocities 
is  the  following: 

.n  .  n  .n  Vp  dvo  .n  dw^  .  q 

61  =  — - hw  Ki,  62  — - - - h  w  1\2,  ^12  —  U  Ki, 

as  Tq  as  as 

.  90  ■  0  dro 

M  ^'2  —  7^- 

as  Tq  as 

The  elastic  deformations  tensor  components  61,62,612  and  the  parameters  A+,  A_  are 
also  assumed  to  be  distributed  along  the  shell  thickness  like  the  components  of  the  full 
deformations  tensor.  One  obtains: 

61  =  6°  -  ZX'i,  62  =  -  2X2,  eu  =  612, 

A+  =  A“  -  2xJ,  A_  =  A°  -  2X^. 

The  tensions  and  moments  are  expressed  in  terms  of  elastic  deformations  components, 
enthropy  and  functions  of  damadge  parameters. 


AO 

AO 

2e°  - 

1  ^ 

1  —  UJ 

1-u^  2e0 

AO 

AC 

2e°- 

UJ 

1  -tu 

1  ^ 
1  ~  UJ 

2eS 

h  —  02361  +  02262  —  d2{r}  —  rjo)  — 


AC 

g  =  20666°,  -  3-^0^, 

Here  the  following  notations  are  used: 

e?,  =  \/wr  +  (e§r-<!M+3«r.  m = b+«a»  . 


a,iij  are  damadge  parameters  and  0,  A,  A,  O  are  constants. 


The  damadge  parameters,  having  zero  values  at  initial  state  of  undamadge  material, 
are  evaluated  using  the  following  differential  kinetic  equations. 


a 


\  —  UJ 
,0 


1  —  u 


oO 


l-uo  /  Vl-o; 


—  e. 


\  /  A° 

=  - aAhI  - A*  1  . 


1-uj  7  \^i  -tj 

Here  H{x)  is  the  Heavyside  function,  T),e*,e^,A*  are  constants  and 


A°  =  7(A0)2  +  (A0)2. 

The  elastic  deformations  tensor  terms  and  terms  of  A+,  A_  which  were  defined  above, 
are  evaluated  by  the  following  set  of  kinetic  differential  equations. 


gO  ^  +  F 


-  -  /?2  9-  +  -  m)\ 


j^ACa  el  +  el  _ 
1  —  a;  1  —  a;  2e° 

AQiio  AC(x  T  6^ 


e°  =  4  -  -  C^e^2  +  Ff^  +  F^^^^  -  +  d.(v  -  ^); 


— Xi  —  —ki  +  C'^Xi  +  C2X2  —  F 


y4C'a'Xi+X2  V  o-  V 


-X2  —  —^2  +  C2X1  +  F'iX2  ~  F- 


I-UJ  2el 

ACa  xi  +  X2 


UJ 


2e?, 


p2  ql  -  A  ?-  +  daVx\ 

P2(l-Fd^fi^\ 


ACn 

e?,  =  e?,  -  2Q'66C'66eL  +  3a, 


A°  —  —Cl  (e^  +  63)  +  2^2  - T  "*■  P2 


KQu 


66  Q  ) 

l-u  el 

ACa  6^  +  62  0 


1  -  a;  '  1  -  a;  26° 


+  79+  -^0); 


. .  - 


-X+  —  Fi  (xi  +  X2)  -  p2  Y 


A"  --C27e^-e")  +  3/?2^ 


P  „0', 


u  2e0 

ACa  el  —  63 
1  —  u  2e?. 


+  379  ; 


■-A  \  o  /3+  ACa  Xi  X2  r>.,„x 


The  expressions  above  use  the  following  notations. 


-379  A 


Constants: 

Cl  =  ai2C'i2  +  0'22C22  +  0^23^23,  C*^  =  ai2C’i2  +  a23C22  +  0:22^23, 
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Variables: 


Cl  —  Pi  (7i2  +  p2  (^22  +  C'23),  C2  —  C22  C231 

F  =  a'22  +  Q-'23i  da  =  ai2di  +  Fd2,  djs  =  Pidi  +  2/?^ d2. 


ql=.B^{Al-9{T-To))-^^ 


A/^Dua  a®  0  r.  A  0  Aa^<^a  A° 


a;  AO 


gO  =  B_A^_  - 


AaDua  xi 


V  r>  /  A\  I  A4-  Y  D  /  A\  , 

=  B+{-x+)  +  - - 9-  =  B^{-X-)  + 


\~L0  AO’ 
AaDua  X- 


l-o;  AO^ 

,  0  _  D  ,A 


\-LO  AO’ 

^  =  5+0AO  ,  rj^  =  B+0X+- 

The  enthropy  r;  averaged  by  the  thickness  of  the  shell  is  evaluated  from  the  following 
eciuation. 

rjT  =  Ni(e^  —  6°)  +  ^2(^2  ~  ^2)  T  ~  ®i«)  ~  A'j.g^  —  A_g°  +  Acj^  +  Ad^  +  AaA^. 

The  temperature  is  evaluated  by  elastic  deformations  and  enthropy: 

T  —  Tq  —  d2{e\  +  eg)  +  P{r]  —  rjo). 

If  we  regard  axisymmetrical  closed  shell,  then  there  are  at  least  2  places  of  the  shell 
laying  on  the  axis.  For  the  points  belonging  to  the  axis  some  expressions  and  equations 
change. 

The  expressions  for  curvatures  Ki  and  /Fg  will  be: 

d^XQ  drp 

Kg  =  AT  - 


drci 

9s 


The  equations  of  motion  change  to: 


u  ~  0 


cFitP  3C)  1 

=  2^  -  iViATi  -  K2K2  + -An(-A/2) 

aF  as  h 

0-0 

Linkage  between  the  velocities  of  deformations  and  velocities  of  displacements  changes 


to: 


diP 


dip 


e?  =  ^  +  w^Ki,  =  —  +  w^K2,  4  =  0, 


k  ~k 

Ai,l  —  /Cg  —  . 

as 
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One  can  note  that  the  following  linkages  must  take  place  on  the  axis: 

■^1  =  -^2)  Q  =  0)  Ml  =  M2,  Cl  =  €2,  €iz  =0,  fci  =  k2- 

The  present  mathematical  model  enables  to  determine  the  parameters  of  an  axisym- 
metrical  composite  damageable  shell  under  the  influence  of  internal  loading. 


1.3  Numerical  modelling. 


The  mathematical  model  described  above  was  used  to  construct  numerical  model. 
This  numerical  model  is  based  on  finite  elements  approach,  but  since  the  equations  contain 
time-derivative  terms  in  significantely  complex  nonlinear  form,  the  classical  finite  elements 
approaches  using  variational  principles  cannot  be  feasibly  applied  to  the  problem.  We  used 
the  so-called  semi-discrete  Galiorkine’s  method  [7]  with  some  modifications.  Details  of 
calculations  are  listed  below. 

We  have  got  the  following  vector  of  variables,  changing  with  time  and  space; 

W  —  (h°,m°,©,u°,m°,©,e?,e2,e°^,  A;i,A:2, 

e?,  4, 4z,Xi ,  X2,  A°  ,  ,  x^,  xj ,  Oi,  cua,  r?). 

Each  member  of  the  vector  W  has  a  governing  differential  equation  evaluating  it  with 
time,  so  the  whole  set  of  differential  equations  can  be  expressed  as: 

dW 

^  =  dW{W). 

The  set  of  expressions  dW{W)  uses  some  significant  variables  not  beloning  to  W  (no 
special  differential  equations  for  them).  We  can  collect  them  in  the  vector  V\ 

K  =  V(W,P)  =  {x'’,r\KuI<2,N^,N2,Q,Mi,M2,T). 

Note  that  P,  i.e.  loading  on  the  wall,  is  obtained  as  an  external  parameter  or  taken  from 
some  external  calculations  (hydrocode  application  for  example),  and  it  is  used  to  calculate 
the  vector  V. 

The  vector  W  consists  of  2  sub- vectors:  W  =  (Wi,  1^2)  so  that  the  first  sub-vector  Wi 
contains  variables,  which  governing  equations  contain  space  derivatives  in  the  right  side. 
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The  equations  for  W2  variables  do  not  contain  such  derivatives  explicitly.  The  sub- vector 
Wi  is  the  following: 

ITi  =  (h°,w°,0,e?,e2,e?^,fci,A:2). 

The  finite  elements  technique  is  used  to  evaluate  expressions  for  dWi. 

The  axisymmetrical  shell  is  represented  by  a  computational  domain  with  Lagrangian 
coordinates  (s,  z).  In  our  case,  variables  change  only  with  the  coordinate  s.  As  we  do  not 
regard  branches  in  this  approach,  then  the  domain  is  a  segment  sm[0,5'].  The  value  of 
S  is  the  length  of  the  generating  line  from  one  pole  of  the  shell  to  another.  We  cut  the 
domain  into  elements,  which  are  segments  [sj,Sj+i],  where  0  <  i  —  2,  so  obtaining 

N  —  I  elements  bounded  by  N  nodes. 

Each  variable  is  assumed  to  have  linear  distribution  along  each  element,  so  that  the 
internal  values  are  expressed  linearely  by  the  values  in  the  nodes: 

v{s)  =  Uj  H - (s  -  Si), 

Si+l  -  Si 

where  v  is  any  one  of  the  variables,  s  G  [Sj,  Sj+i],  and  Uj  is  value  in  the  Tth  node. 

The  equations  evaluating  vector  Wi  can  be  represented  in  more  detailed  form: 

dWi  _Edr»  aF 

dt  ds  ds 

We  regard  any  node  i  and  the  elements  adjacent  to  it.  In  our  case,  there  are  generally 
2  elements  placed  to  the  left  and  to  the  right  except  for  the  leftmost  and  rightmost  nodes 
both  laying  on  the  axis.  We  integrate  the  equations  on  halves  of  each  adjacent  elements, 
and  using  our  assumption  of  linear  distribution,  each  integral  can  be  evaluated  analytically 
thus  simplifying  calculations.  As  there  are  no  second-order  derivatives  in  the  right  side 
of  the  equations,  then  the  assumption  of  linear  distribution  is  enough  to  obtain  correct 
expressions.  For  the  right  side  we  obtain  an  expression  depending  on  3  successive  nodal 
values  of  W  and  V  as  result  of  integrating.  Integrating  the  left  side  and  equating  it  to 
right,  we  obtain  implicit  set  of  linear  ecpiations  with  the  same  matrix  evaluating  dWi.  In 
our  case  matrices  of  all  the  eciuations  are  3-diagonal,  and  this  helps  us  to  invert  them 
with  direct  method  of  binary  substitution. 

The  complete  operation-of  transfer  from  the  time  layer  n  to  the  next  layer  consists  of 
following  steps. 

a.  Determine  the  timestep  value  ht  according  to  the  Courant  criterion  modified  to 
avoid  overshoots  in  evaluating  kinetic  equations. 

b.  Get  the  values  of  loading  at  the  current  time  moment. 

c.  Evaluate  E”  =  E(1T",P"). 
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d.  Evaluate  dW^  consisting  of  dVE”  and  dW^  ■  For  the  first  set  of  equations  we  apply 
finite  elements  techniciue,  the  second  set  is  simply  evaluated  from  expressions. 

e.  Predictor  -  preliminary  evaluation  of  W  to  the  next  time  layer; 

W  =  VF”  +  ht  ■  dW^. 

f.  Increase  time  counter  and  get  the  loading  on  the  next  layer. 

g.  Evaluate  V  =  F(VF,F"). 

h.  Evaluate  dW  using  the  same  technicpie  as  in  the  item  d. 

i.  Corrector  stage  -  obtain  values  of  W  on  the  next  time  layer: 

Vl/n  +  l  ^ 

j.  Filtering  the  probable  non-physical  oscillations  of  velocities  of  displacements.  This 
stage  was  suggested  in  [7]. 
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Chapter  2 

Part  II.  Results  of  numericcJ 
modelling. 


2.1  Foreword. 


The  present  Part  contains  the  description  of  the  results  of  numerical  modelling  of  dyna¬ 
mical  deforming  of  the  walls  of  an  axisymmetrical  containment  made  of  thermoviscoelastic 
laminated  damageable  composite  material  under  the  influence  of  internal  loading.  Two 
types  of  internal  loading  were  regarded  to  demonstrate  the  utility  of  the  composite  model. 
The  first  one  was  characterised  by  a  uniform  increase  of  pressure  in  the  containment  that 
brought  to  an  expansion  of  the  shell  and  accumulation  of  damages  wherein  nonuniformity 
was  introduced  only  by  the  initial  difference  of  curvatures  in  different  zones  of  the  shell.  In 
fact  it  was  shown  that  for  the  case  of  uniform  internal  loading  of  a  shell  having  the  shape 
different  from  a  spherical  one  there  appeared  nonuniformities  in  accumulation  of  damages 
and  specific  dissipation  was  maximal  in  the  zones  of  maximal  curvature  gradients. 

The  second  type  of  loading  regarded  was  the  one  in  reflected  shock  waves  caused  by 
internal  explosion  in  a  gas- filled  containment.  The  rates  of  wall  loading  were  determined 
having  incorporated  the  present  composite  shell  model  into  a  hydrocode.  To  garantee  the 
axisymmetrical  loading  pattern  the  centre  of  the  explosion  could  be  placed  in  any  place  of 
the  axis  of  symmetry  of  the  shell.  A  more  detailed  description  of  the  model  determining 
the  wall  loading  and  of  the  flow  pattern  in  the  gas  inside  the  containment  will  be  given 
below. 
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2.2  Determining  the  internal  loading  of  a  contain¬ 
ment. 


The  internal  nonuniform  loading  of  the  walls  of  the  containment  was  determined  within 
the  frames  of  an  assumption  that  the  containment  was  filled  in  with  an  inert  or  chemically 
reacting  gas  mixture.  The  internal  explosion  was  simulated  by  introducing  discontinuity 
in  the  initial  conditions:  assigning  high  pressure  and  density  values  for  the  parameters 
inside  a  relatively  small  zone  on  the  axis  of  symmetry.  To  describe  the  originating  flow  a 
gasdynamical  model  for  turbulent  flows  in  confined  volumes  was  used. 


2.2.1  Mass  balance  equations  in  the  gas  phase. 


The  mass  balance  of  the  /c-th  component  is: 

dt{pk)  +  V  •  {pkUk)  =  Mk, 

where  pk  is  the  density  of  the  component  (mass  per  gas  phase  volume  unit),  Uk  is  the 
velocity  vector  of  this  component  and  .Adj.  is  mass  flux  to  the  /c-th  component  from  the 
other  components. 

The  following  notations  and  definitions  are  introduced:  p  =  J2pk  for  gas  phase  density, 

k 

pu  =  pkUk  for  gas  phase  general  velocity  vector,  Wk  —  Uk  —  u  for  the  velocity  of  /c-th 

k 

component  diffusion,  h  =  PkWk  for  the  diffusive  flux  of  A;-th  component,  Yk  =  pk/p  for 
mass  fraction  of  the  A>th  component  so  that: 

=  E4  =  o. 

k  k 

With  these  definitions  the  mass  balance  equation  for  the  fc-th  component  in  the  gas 
phase  can  be  transformed  into: 

a,{fiY,)  +  V  ■  (pci-i)  =  M,-v-h.  (1) 
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After  summing  up  the  equations  (1)  for  all  the  components  the  general  mass  balance 
equation  for  the  gas  phase  is  obtained: 

dt(p)  +  V  ■  (pm)  =  0.  (2) 

After  Favre  averaging  of  the  ecjuation  (2)  with  the  p  weight  the  averaged  mass  balance 
equation  for  the  gas  phase  is  obtained: 

dt{p)  +  V  •  {-pu)  =  0.  (3) 

By  averaging  (1)  in  the  same  way  and  using  the  notation  II  =  pu'^Yk”  for  the  turbulent 
diffusion  flux  the  averaged  mass  balance  for  the  k-th  component  in  the  gas  phase  is 
obtained:  _ 

dtipYk)  +  V  •  (puYk)  =  -  V  ■  (4  +  4).  (4) 


2.2.2  Momentum  equation  for  the  gas  phase. 

For  the  k-th  component  in  the  gas  phase  the  momentum  equation  is: 

dt[pkUk)  +  V  •  {pkUk  0  Ilk)  =  —'^Pk  +  PkQ  +  V  •  Tfc  +  ICk, 

where  Pk  is  the  partial  pressure  in  the  component,  r*.  is  the  partial  viscous  tensor,  g  is 
the  gravity  acceleration  vector  and  )Ck  is  the  momentum  flux  to  the  component  from  the 
other  components. 

Summing  up  these  equations  for  all  components,  the  following  equation  is  obtained: 

dt{pu)  +  V  •  {pu  <S)u)  =  - Vp  +  +  V  •  r,  (5) 

where:  p  =  J^Pk  is  the  pressure  in  the  gas  phase  according  to  Dalton’s  law, 

k 

k 

is  the  condition  of  agreement  for  momentum  fluxes  from  all  the  components.  It  is  also 
assumed  that  the  viscous  stresses  tensor  for  multicomponents  flow  is  r  =  Y^i^k^  ^Ik®h)- 
Following  the  Favre’s  techniques  [4]  the  averaging  of  (5)  yields: 

dt{pu)  +  V  •  {pit  u)  =  pg  —  Vp  +  V  •  ("^  +  (6) 

where  r*  =  —pu”  ®  u"  is  the  Reynolds  stress  tensor. 
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2.2.3  Energy  balance  equation  for  the  gas  phase. 


For  each  component  in  the  gas  phase  it  is  possible  to  write  the  following  energy  balance 
equation: 

dt{pkEk)  +  V  •  [pkUkEk)  =  PkUk  -g-V-  (pkUk)  -  V  •  V  +  V  •  (r^  •  +  pkOk  +  ^k, 

where  Ek  —  Ck  +  \'Ei  is  a  sum  of  internal  and  kinematic  energy  for  the  A;-th  component, 
Iqk  is  conductive  heat  flux  to  the  k-th.  component,  Ok  is  the  radiation  heat  flux,  £k  is  the 
heat  flux  from  the  other  components.  It  is  assumed  that  the  gas  phase  is  a  mixture  of 
perfect  gases,  so  that; 

—  CykT  hokj  Pk  PkE  1 

where  T  is  the  temperature,  Cyk  is  the  heat  capacity  at  constant  volume,  Wk  is  the  molar 
mass  and  Rg  is  the  gas  constant, /lo*,  is  the  chemical  energy  of  the  fc-th  component. 

The  following  notations  for  the  mixture  energy  and  generalized  heat  flux  are  used: 

£  =  ^n£.  =  E(n-e'.  +  ^)  +  ^. 

k  k  ^Pk 

4  =  YLV*  +  -  (n  -  4l  =  4  +  E  '*‘4. 

k  ^Pk  Pk  k 

where  ^  is  the  enthalpy  of  k-th  component.  With  these  notations  the  summing 

up  of  the  energy  balance  equations  for  all  components  leads  to  the  following  equation: 

dt{pE)  +  V  •  (puE)  =^pu-g-V-  (pu)  -V  ■  Iq  +  V  ■{t-u)+  pO,  (7) 

where  the  condition  of  agreement  for  energy  exchange  between  the  components 

E«^  =  o 

k 

is  taken  into  account. 

The  terms  of  the  second  and  higher  orders  of  diffusion  fluxes  are  neglected  in  the 
expression  for  the  mixture  energy.  The  equation  of  state  for  the  gas  phase  and  the 
expressions  for  the  internal  energy  and  enthalpy  then  take  the  form: 

p  =  B,prEi^.  e  =  EU(c„»r  +  /i„t).  ft  =  Eu(cpi.r  +  ft„o.  (s) 

k  ^  k  k  k 

Averaging  (7)  by  Favre’s  technicpies  allows  to  obtain: 


Assuming  the  expression  (8)  to  be  true  also  for  the  averaged  values,  brings  the  following 
assumptions: 

E  E  c.ier'U” « e,  E 

k  P  k  k 

With  these  assumptions  finally  the  averaged  energy  balance  equation  can  be  obtained: 

dtipE)  +  V  ■  (puE)  =-pu-  g-V  ■  {pu)  -  V  •  (7,  +  /J)  +  V  •  (t^  +  r*  •  m)  +  p9.  (9) 

The  notation  E  =  pJ2cpkYku”T'  +  Yl{CpkT  +  hok)Ik  “  +  I7(CpfcT  +  hok)Ik  is  used 

k  k  .  ^  . 

for  the  generalized  turbulent  heat  flux.  The  expression  for  the  mixture  energy  E  = 
E  YkiCvkYEhok)  contains  the  averaged  kinematic  energy  of  turbulent  pulsations 

fc(the  last  term  is  usually  neglected  in  calculations  because  it  is  assumed  to  be  much 
smaller  than  the  scalar  product  of  the  averaged  gas  velocity) . 


2.2.4  Fluxes  modelling  and  closing  the  K-epsilon  model. 


According  to  the  standard  A:-epsilon  model  for  compressible  gas  flows,  the  following 
model  for  turbulent  fluxes  is  introduced: 

2  2 

T  eE  =  pY  E  E){Vu  E  VtE  -  -V  •  ul)  -  -pkl, 

o  o 

hEll^-p{DE-)V-Y, 

YeJ^^  -  E  +  -)V  •  T,  (10) 

k 

where  I  is  the  unit  tensor  of  the  second  order.  The  energy  dissipation  term  in  the  equation 
(9)  can  be  transformed  as: 

V  •  E  E  •  u)  —  Y  ■  ((r  +  E)  ■  u). 


The  turbulent  kinematic  viscosity  E  is  modelled  according  to  A:-epsilon  model  as: 


E  =  C, 


E 


(11) 


The  model  is  closed  then  by  2  equations  for  the  kinematic  energy  of  turbulent  pulsations 
k  and  its  decay  due  to  dissipation  e: 


dtipk)  E  V  •  {pPk)  =  V  •  {pY  E  —)Vk)  eE  :Vu-  pe, 

O'k 


(12a) 
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(126) 


a,(p£)  +  V  ^  (Tmi)  =  V  ■  {p(v  +  -)Ve)  +  ■.  Vi  -  C^pc). 

(jf  k. 

In  this  equations  the  Reynolds  tensor  r*  can  be  modelled  as; 

=  -pu\Vu  +  ■  ul)  -  Ipkl. 

U  O 

The  term  containing  k  in  the  last  expression  is  negligible  as  well  as  in  the  expression  for 
the  sum  of  viscous  and  Reynolds  tensors. 

The  constants  take  the  following  values; 

C;  =  0.09,  (j,.  =  1,  ae  =  1.3,  Cu  =  1.45,  C'ae  =  1.92,  =  1.  (13) 

The  turbulent  Prandtl  and  Schmidt  numbers  are  assumed  to  be  equal  to  1  in  (13). 


2.2.5  Initial  and  boundary  conditions. 


The  initial  mean  flow  velocity  in  the  air  was  supposed  to  be  zero  u  =  0,  but  the  initial 
level  of  turbulence  was  introduced;  k[t  =  0)  —  ko]e{t  =  0)  =  Cq.  A  definite  volume  on  the 

axis  of  symmetry  (0  <  ?’  <  7’o)  ,  Xi  <  x  <  X2  was  considered  to  be  occupied  by  dense 

reaction  products  under  high  pressure.  The  rest  of  the  containment  was  considered  to  be 
occupied  by  the  gas  of  ambient  pressure  and  density.  At  t  =  0  the  motion  starts  and  the 
initial  discontinuity  in  the  pressure,  density  and  temperature  fields  brings  to  formation  of 
shock  waves  propagating  inside  the  containment  and  being  reflected  by  the  walls. 

The  boundary  conditions  on  the  walls  of  the  containment  are  the  following; 

u  =  0  —  the  no-slip  condition; 

dT  dYk  dk  de  ^  ^ 

-v—  =  -v —  =  —  =  —  =  0  —  the  Neumann  conditions. 

UZ  UZ  OZ  UZ 


2.2.6  Computational  techniques. 


The  calculations  are  carried  out  using  two  demi-steps  to  determine  new  parameters 
for  the  gaseous  phase. 

The  value  of  the  timestep  is  recalculated  on  the  base  of  Courant  criterion. 
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The  system  of  gasdynamics  equations  can  be  rewritten  in  the  following  vector  form; 

(14) 


where 


U  - 


d  d  ^ 

d 

-^U  F  -^E  +  ^F+  =  H, 

at  ax 

or 

"  1 

^  pu  ^ 

f  0  ^ 

pu 

pu^  +  p 

'^xx 

pv 

puv 

Tcz 

pE 

,  E^E^  +  E^  = 

puE  +  up 

+ 

'^'^xx  '^'^xz  T  Jqx 

pk 

puk 

Ex 

pe 

pu€ 

Ex 

V  p^'A-  y 

\  puYk  y 

^  Ex  / 

^  pv  ^ 

^  0  ^ 

puv 

"Pcz 

pv"^  +  p 

pvE  +  vp 

-1- 

UT^;^  -f- 

pvk 

Ez 

pvt 

Ez 

\  pvYk  i 

\  Ez  / 

H  =  H^  = 


Pk 


-pe 
liCuPk  -  C2epe) 

Mk 


The  formulae  above  contain  the  following  notations:  x^r  -  axial  and  radial  coordinates  re¬ 
spectively,  u,  V  -  mean  velocity  components  in  this  coordinates,  Pk  is  the  turbulent  energy 
production  term  (equation  (12)).  The  superscript  ” H”  means  convection  and  production 
terms  (hyperbolic  part),  ”F”  -  generalized  diffusive  terms  (including  viscous,  diffusive  and 
thermoconductive  terms,  (parabolic  part)).  Note  that  the  system  (14)  is  mathematically 
overdefined,  for  ^'1*  =  1)  but  in  numerical  calculations  one  more  equation  is  kept  for 
the  control  of  precision  and  correction  of  results. 

To  solve  the  system  (14)  splitting  by  coordinates  is  used;  it  simplifies  the  solution  and 
increases  the  possible  timestep.  The  splitting  is  made  according  to  MacCormack  [5].  This 
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splitting  represents  the  general  operator  L{At)  transferring  t/”  to  the  next  timestep  into 
in  the  form; 

L{At)  =  Lr{At,;)Lr{Atr)LriAtr)Lx{At^), 
or  L{At)  =  Lr{Atr)L^{At^)L^{At^)LriAtr).  (15) 

The  sequence  of  operators  in  (15)  yields  the  condition  of  symmetry.  To  yield  the  condition 
of  timesteps,  which  together  with  the  condition  of  symmetry  ensures  the  second  order  of 
approximation,  it  is  necessary  to  have: 

At  =  2A4  =  2Atr, 

so  that  Atx  =  Atr  =  At/2.  Both  sequences  of  operators  in  (28)  are  able  to  represent 
the  general  operator  L{At)]  to  avoid  the  accumulation  of  disagreements  the  sequences  are 
changed  at  each  timestep. 

Each  operator  itself,  and  Lr,  is  also  split  into  two  parts:  one  part  is  parabolic, 
or  L^,  the  other  is  hyperbolico-parabolic,  or  ,  and  contains  terms,  which  are 
not  present  in  or  .  This  gives  us  the  MacCormack  rapid  solver  techniques  [5]. 
The  parabolic  operator  is  solved  using  the  implicit  Laasonen  scheme,  the  hyperbolico- 
parabolic  operator  is  solved  using  FCT  techniques  [6].  This  splitting  removes  viscosity 
from  the  timestep  criterion  and  reduces  it  to  the  Courant  criterion. 


2.2.7  Results  of  numerical  modelling  of  internal  loading. 


The  results  of  numerical  modelling  of  shock  waves  propagation  and  reflection  inside 
the  containment  are  shown  in  Figs.  2-3.  The  figures  show  the  density  and  pressure  fields 
inside  the  containment  for  different  times.  The  length  of  the  cylinder  was  Im,  the  radius 
-  0.5m.  The  map  of  characteristic  values  of  parameters  is  given  in  the  upper  right  parts 
of  the  figures.  The  characteristic  times  are  also  given  on  the  right  screen.  The  black  line 
segments  on  the  figures  reflect  gas  velocities. 

The  initial  volume  where  the  explosion  takes  place  can  be  seen  in  Fig. 2. a  illustrating 
the  initial  density  distribution  at  f  =  0.  Expansion  of  hot  reaction  products  brings  to 
a  formation  of  a  ball-shape  zone  of  a  compressed  gas  (Fig. 2. 5.)  that  after  the  reflection 
from  a  lower  wall  comes  to  a  hemispherical  shape  (Figs.2.c,  d.).  The  converging  rarefac¬ 
tion  waves  bring  to  a  decrease  of  density  behind  the  diverging  hemispherical  shock  wave 
(Figs. 2. e,  /).  On  reflecting  the  shock  from  the  side  walls  of  the  containment  density  incre¬ 
ases  behind  the  reflected  waves  {Figs.2.g,h).  The  position  of  the  reflected  waves  can  be 
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easily  tracked  by  the  gradients  of  the  velocity  field.  The  reflected  shock  waves  turn  to  be 
steeper  and  overtake  the  leading  shock  near  the  walls  (Figs.2. h,i^k.)  that  finally  brings 
to  a  formation  of  a  Mach  stem  (Figs. 2. 1,  m).  The  reflected  shock  waves  are  converging  to 
the  axis  and  bring  to  an  increase  of  density  after  their  collision  (Figs.2.A;, /,  m). 

The  zone  of  high  density  originates  near  the  upper  wall  after  the  reflection  of  the 
leading  shock  and  the  secondary  shocks  (Figs.2.?n,  o).  The  reflected  shock  wave  pro¬ 

pagates  faster  near  the  walls  than  in  the  center  (Figs.2.p,  s.)  because  it  meets  in  the 
center  a  more  dense  gas  with  a  large  accumulated  axial  velocity. 

The  pressure  field  evolution  shown  in  Figs. 3. a  —  n  gives  an  idea  of  the  internal  loading 
of  the  shell  of  the  containment.  The  ball-shape  shock  wave  is  formed  at  the  very  beginning 
(Fig. 3. a).  In  reflection  of  the  shock  wave  from  the  nearest  lower  wall  the  loading  starts 
that  spreads  from  the  axis  of  symmetry  (Figs.3.6,  c,  d,  e).  The  converging  rarefaction 
waves  gradually  decrease  the  pressure  in  the  center  (Figs.3.d,  e.)  and  a  high  pressure 
zone  is  present  only  in  the  primary  shock  wave  (Fig.3.e).  The  intensity  of  the  primary 
shock  wave  decays  but  a  new  pressure  increase  takes  place  in  the  transverse  shock  waves 
reflected  from  the  side  walls  (Figs. 3./, 5').  Converging  shock  waves  bring  to  a  pressure 
increase  on  the  axis  of  symmetry  (Figs.3./i,  i).  On  reflecting  the  primary  shock  from  the 
upper  wall  the  pressure  increases  in  the  reflected  shock  wave  (Figs.3.A;, /,  m.)  and  then 
gradually  decreases  due  to  the  influence  of  rarefaction  waves  (Fig.3.n). 

This  typical  example  of  wall  loading  in  internal  explosion  gives  wide  possibilities  to 
investigate  the  behaviour  of  the  model  for  damageable  composite  shell  in  nonuniform 
loading. 
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2.3 


Dynamical  deforming  of  walls  and  accumulation 
of  damages. 


Here  we  describe  the  results  of  model  calculations  of  shell  parameters  under  the  in¬ 
fluence  of  internal  loading.  The  characteristic  values  of  material  parameters  are  given  in 
the  Table  1.  The  two  cases  of  loading  -  uniform  and  nonuniform  -  are  regarded.  The 
characteristic  value  of  the  critical  dissipation  was  taken  arbitrary.  The  model  incorporates 
as  well  three  additional  material  constants  characterizing  the  resistance  of  the  laminate 
to  accumulating  the  damages  in  shear,  tension  and  delamination.  The  values  of  those 
parameters  were  also  chosen  arbitrary  due  to  the  lack  of  kinetic  data  on  accumulation  of 
damages  in  the  laminates. 


2.3.1  Evolution  of  the  shell  parameters  in  uniform  dynamical 
loading. 

The  simplified  case  of  uniform  internal  loading  was  regarded  to  investigate  the  expansion 
of  the  shell  and  accumulation  of  damages  wherein  nonuniformity  was  introduced  only  by 
the  initial  difference  of  curvatures  in  different  zones  of  the  shell.  The  internal  loading  was 
organized  due  to  uniform  pressure  increase  inside  the  containment  up  to  Pmax  —  lO^Pa 
with  the  characteristic  time  10“®s.  The  initial  thickness  of  shell  was  h  —  initial 

radius  -  1???.,  initial  length  -  2m. 

Fig.-^  shows  the  initial  shape  of  the  containment  in  Euler  coordinate  system  and 
Fig. 5  -  in  a  Lagrangian  one.  The  characteristic  internal  loading  profile  for  the  present 
case  is  shown  in  Fig. 6.  The  shape  of  the  deformed  containment  for  the  characteristic  time 
5.95  ■  10~^s  is  shown  in  Fig. 7.  It  can  be  seen  from  the  Fig. 7  that  the  shell  expanded  and 
changed  the  curvatures. 

Figs. 8  a  —  b  show  the  curvature  of  the  generating  line  of  the  cylinder  for  two  successive 
times.  Figs. 9  a-b  show  the  curvature  of  the  orthogonal  geodesical  line  for  the  same  times. 
It  can  be  seen  from  the  figures  that  both  curvatures  have  the  tendency  to  change  (increase 
or  decrease)  towards  the  mean  value. 

Figs.  10  a—b  show  the  displacements  of  the  shells  elements  in  the  directions  orthogonal 
to  the  surface  (along  the  lagrangian  axis  z).  The  vertical  lines  in  the  Figs.  10  and  successive 
figures  indicate  the  places  of  initial  discontinuities  of  curvatures.  The  local  maxima  of 
displacements  w  are  found  in  the  zones  of  the  minimal  initial  curvature  that  corresponds 
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with  the  tendency  of  the  containment  to  come  to  a  spherical  shape  under  the  influence  of 
a  uniform  loading. 

Figs. 11  a  —  b  show  the  displacements  of  the  shell’s  elements  u  along  the  generating 
line  in  lagrangian  coordinates  s.  The  displacements  along  the  s-axis  have  different  signs 
for  different  elements.  The  positive  or  negative  displacements  are  determined  by  the 
initial  positions  of  the  elements  in  respect  to  the  boundaries  and  the  middle  point.  The 
displacements  of  the  elements  located  on  the  axis  of  symmetry  (boundaries)  and  in  the 
middle  zone  are  ecjual  to  zero.  The  inverse  symmetry  in  respect  to  the  middle  point  is 
also  evident.  These  results  qualitatively  testify  the  accuracy  of  the  worked  out  model. 

Figs.  12  a—i  illustrate  the  velocity  profiles  for  the  orthogonal  displacements  for  different 
times.  It  is  seen  that  at  the  initial  stage  of  the  process  velocities  of  displacements  change 
their  values  and  signs  very  rapidly.  The  behaviour  of  the  velocities  shows  the  oscillating 
character  of  the  expansion  of  the  shell  in  rapid  dynamical  loading.  Later  the  attenuation 
of  the  oscillations  takes  place  due  to  viscous  dissipation  and  the  expansion  velocity  is 
stabilized. 

Figs. 13  a  —  m  illustrate  the  velocity  profiles  for  tangential  displacements  along  the 
generating  line  for  different  times.  The  results  of  numerical  modelling  evidently  illustrate 
the  tangential  oscillations  of  the  elements  and  nonuniformities  in  velocity  profiles  caused 
by  the  initial  difference  of  curvatures  in  different  zones  of  the  shell.  Thus  even  in  a  uniform 
internal  dynamical  loading  of  shells  with  variable  curvatures  there  appear  nonuniformities 
in  the  parameters  for  different  elements  of  the  shell. 

Figs. 14-16  show  the  profiles  of  tangential  strains  in  both  orthogonal  geodesical  direc¬ 
tions  of  the  shell  6°  and  the  strain  for  one  and  the  same  time.  It  is  seen  from  Fig. 14 
that  Cl  reaches  its  maximal  value  in  the  zones  of  largest  curvature  of  the  generating  line. 
The  Fig. 15  shows  that  the  maximal  values  of  can  be  found  in  the  middle  zone  wherein 
the  radial  displacements  are  also  maximal.  And  local  extrema  of  the  (Fig.  16)  are  in 
the  zones  of  initial  discontinuity  of  the  curvature  of  the  generating  line  of  the  shell.  The 
above  results  correspond  well  with  the  physical  nature  of  the  phenomenon. 

Figs.  17-19  show  the  profiles  of  the  corresponding  elastic  strains  e°,e2,  for  the  same 
times.  It  can  be  seen  that  by  that  time  the  difference  between  the  elastic  and  total 
deformations  has  already  become  significant. 

Figs. 20-21  show  the  profiles  of  parameters  A+  and  A_  characterising  the  difference  of 
composite  components’  strains  in  the  layers  plane  [2]. 

Figs. 22-24  show  the  profiles  of  damage  parameters  uj,a,uj^  reflecting  accumulation  of 
damages  in  tension,  shear  and  delamination.  The  zones  of  maxima  of  damages  in  tension 
coinside  with  that  for  delamination  while  the  zones  of  maxima  of  damages  in  shear  are 
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located  near  the  places  of  initial  discontinuity  of  curvature  of  the  generating  line  of  the 
shell. 

Figs. 25  a  —  b  illustrate  the  Ni  force  in  the  shell  along  the  generating  line  for  the 
two  successive  times.  Figs. 26  a  —  b  illustrate  the  tangential  to  the  shell  N2  force  in  the 
orthogonal  direction  for  the  same  times.  Figs. 27  a  —  c  show  the  profiles  of  the  force  Q 
orthogonal  to  the  tangent  plane. 

Fig. 28  illustrates  the  temperature  increase  in  the  shell  due  to  irreversible  processes.  At 
the  initial  instant  it  has  the  zero  level.  It  is  seen  that  the  maximal  increase  of  temperature 
can  be  found  in  the  zones  of  maximal  gradients  of  curvature  of  the  generating  line. 

Figs. 29  a  —  d  show  the  accumulation  of  the  specific  dissipation  D  in  the  shell  for 
successive  times.  The  growth  of  specific  dissipation  is  maximal  in  the  zones  of  maximal 
gradients  of  curvature  of  the  generating  line  and  in  the  zones  of  maximal  deformations  as 
well  (the  midle  point). 

The  results  show  that  for  axisymmetrical  shells  under  the  influence  of  a  uniform  in¬ 
ternal  loading  the  damaged  zones  are  likely  to  appear  in  the  zones  of  maximal  curvature 
gradients  and  in  the  zones  of  maximal  expansion.  Thus  the  presence  of  angles  essentially 
increases  the  probability  of  a  breakup  contrary  to  the  spherical  shape  of  the  containment 
that  is  proved  to  be  optimal  for  the  cases  of  uniform  internal  loading.  This  coinsidence  of 
the  numerical  results  with  the  conclusions  derived  from  the  existing  theoretical  solutions 
makes  us  sure  that  the  worked  out  model  of  the  damageable  thermoviscoelastic  laminated 
composite  material  adequately  describes  the  physical  processes  in  internal  loading. 
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2.3.2  Shell  parameters  in  nonuniform  dynamical  loading. 

To  investigate  the  behaviour  of  a  shell  in  nonuniform  loading  the  dynamics  of  wall  loading 
was  assumed  to  be  that  determined  in  the  section  2.1.  The  characteristic  size  of  the  shell 
was  assumed  the  following;  length  -  Im,  radius  -  0.5m,  initial  wall  thickness  -  0.5-  lO^^m. 
The  size  of  the  shell  corresponds  with  that  of  the  containment  described  in  2.1.  Thus  the 
present  shell  is  twice  smaller  than  the  one  investigated  in  the  previous  section  in  a  uniform 
loading.  To  keep  the  similarity  the  thickness  was  also  reduced  twice  in  comparisson  with 
the  previous  case.  That  agrees  with  the  static  similarity  criterion  but  of  course  can 
introduce  differences  for  a  dinamical  case. 

The  Figs. 30a  —  m  show  the  internal  wall  pressure  distribution  along  the  generating 
line  of  the  containment  (the  lagrangian  coordinate  s)  for  different  times  for  the  case  of 
nonuniform  loading  described  in  the  section  2.1.  It  is  seen  that  in  the  very  beginning  of 
the  process  there  is  a  sharp  pressure  increase  near  the  bottom  (Figs. 30. a  —  h.)  that  is 
accompanied  by  the  increase  of  pressure  near  the  side  walls  in  reflected  waves  (Figs. 30. c  — 
d).  But  the  high  loadings  are  not  durable  and  pressure  decays  very  quichly  due  to  the 
influence  of  the  rarefaction  waves.  After  the  multiple  reflection  from  the  upper  wall 
(the  right  hand  side  of  the  figure)  and  the  pressure  increase(Fig.30./  —  h)  the  loading 
stabilizes  and  small  axial  oscillations  of  pressure  can  be  observed  (Figs. 30. A:  —  m).  The 
mean  pressure  in  the  containment  finally  is  much  less  than  that  in  the  previous  case  for 
the  uniform  loading. 

Figs. 31.0  —  e  illustrate  the  behaviour  of  tangential  displacement  u  in  nonuniform  loa¬ 
ding  for  the  successive  times.  The  displacements  are  oscillating  nonsymmetrically  in  the 
initial  stage  of  the  process  but  finally  they  form  a  steady  profile  similar  to  that  for  the 
case  of  uniform  loading.  Comparisson  of  the  Fig. 31. e  with  the  Fig. 11. 6  testifies  the  fact. 

Fig. 32.0.  —  e  illustrates  the  changes  in  normal  displacements  w  of  the  shell.  It  is  seen 
that  the  sell’s  elements  start  mooving  in  the  left  hand  side  (the  bottom  wall  of  the  shell) 
and  gradually  all  the  shell  is  being  involved  into  motion.  By  the  time  internal  pressure 
comes  to  a  practically  steady  uniform  distribution  the  normal  displacements  form  a  profile 
(Fig. 32. e.)  similar  to  that  for  the  case  of  uniform  loading  described  above  (Fig. 10. 6). 

The  Figs. 33. a  —  h  and  Figs. 34. a  —  /  illustrate  the  velocities  of  shell’s  elements  in  tan¬ 
gential  to  the  generating  line  and  normal  directions  respectively  for  successive  times.  It  is 
seen  that  relatively  large  elastic  oscillations  present  in  the  beginning  are  being  moderated 
in  the  longrun  by  the  increasing  dissipation. 
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The  dynamics  of  the  dissipation  function  for  the  shell  is  shown  in  Figs. 35. a  —  c.  It 
is  seen  that  the  dissipation  function  has  not  reached  its  critical  value  by  the  time  of  the 
pressure  stabilization  inside  the  containment.  That  happened  due  to  the  fact  that  even 
high  but  short  term  loadings  do  not  bring  much  damages  to  the  material  due  to  the 
relaxations  induction.  Thus  the  model  is  able  to  take  into  account  the  relaxation  time  for 
accumulation  of  damages. 

The  successive  stages  of  the  accumulation  of  damages  in  shear  a  are  shown  in  Figs. 36a— 

/. 

The  last  time  moment  reflected  in  the  figures  actually  is  not  the  last  one  for  the 
continuous  destruction  of  the  material  though  the  mean  internal  pressure  is  rather  low 
already.  The  presence  of  large  differences  in  the  total  and  elastic  deformations  (Figs. 37 
and  38)  shows  that  the  relaxation  process  is  not  finished  yet  and  the  dissipation  will  be 
still  increasing  for  some  time.  But  it  is  not  likely  to  surpass  the  criterion. 

Fig. 39  shows  the  distribution  of  temperature  deviations  in  the  shell  from  the  initial 
state.  The  temperature  increase  in  the  right  hand  side  is  due  to  the  growth  of  dissipation. 
The  decrease  of  temperature  in  the  left  hand  side  is  due  to  cooling  of  the  material  in  elastic 
expansion  that  surpasses  the  influence  of  the  dissipation  growth.  The  local  maxima  of 
temperature  in  the  zones  of  discontinuity  of  curvature  are  due  to  the  influence  of  the 
growth  of  dissipation  that  is  maximal  within  these  zones. 

The  present  example  illustrates  the  fact  that  nonunifomities  of  loading  being  stabilized 
during  a  relatively  short  period  of  time  bring  to  the  deformed  picture  similar-  to  the  case 
of  uniform  loading.  The  intense  short  term  loads  bring  less  accumulated  damages  than 
durable  loads  of  smaller  values. 


114 


ding  on  the  uall,  (Pa) 


SO*  a 


H5 


ding  on  the  uall,  CFa) 


ding  on  the  yatll>  (PaJ> 


SOiC^  • 


Z0.  d. 


HS 


oading  on  the  faiall,  (Pa) 


an  the  uall,  (Fa) 


/■/g,,  30,  gi. 

i2i 


30. 

iZA 


uall,  CPa) 


1Z5 


ding  on  the  wall«  (Fa) 


i27 


lacemetit  tangential  imy 


F‘g.3l.  4. 

/28 


isplacement  tangential  U>  (m) 


;  tangential  U,  (m) 


/30 


ispIaceiBisiit  tan^ntial  U,  Cm) 


(VI)  'f)  iiBUiJou  ‘i^uaeiaovidsi 


i  3  2 1  , 


/32 


isplaceiiient  noriuftl  \4,  Cwt) 

10.0+  I  t 


.  32  <  c  , 


lacement  normeil  U,  Cm) 


/36 


ocity  of  taiigeiitial:  displacemgiift  U,  (m/s) 


,  33 ,  a  , 

1d7 


of  taTigemitial  dispiacement  U«  (m/s) 


locity  of  tangential  displacement  (m/s) 


of  tangential  displacement  U«  Cw^s) 


of  tangemtial  displaceBWgnt  U,  Cw/s) 


A  t  3  3  <  0, . 
ni 


of  tangential  displacement  U,  (m/s) 


/42 


int ial  displaceeieifit  U,  (in/s) 


locity  of  tangcmtial  displacement  U,  (m/s) 


locitjj  of  norml  displacement  U,  (ift<’s) 
Ul,  »10"1 


of  noraiii  displaceiaent  Mj  (m/s) 


i'^  & 


Icwr^ity  of  norwal  displacement  Wj  (iw/s) 


<  3  ^  I  , 


i^7 


displacemimt  y«  (m/s) 


Idcity  of  normal  displacement  Uj  (m/s) 


F,W.  37, 


of  norwail  displacewent  y,  (m/s) 


3^. 


150 


issipated  energy ^  (J/kg) 


SSi  oi . 


iffl 


(J/kg) 


F '  3-5” I  ^ , 


issipated  energy «  (J/Kg) 


Fc^.  35-.  a, 


/53 


anadge  parameter  if^l 


/STS' 


i5G 


IS  7 


teo 


deforinatioii  £11  j,  () 


2>  8 , 
i61 


emperature  deviatioiij  (K) 


References. 


1.  N.N. Smirnov,  A. B. Kiselev,  I.D.Dimitrienko,  V.R.Dushin,  V.F. Nikitin.  Computa¬ 
tional  damage  model  for  composite  materials.  Mathematical  Model  of  Dynamical 
Deforming  and  Fracture  of  Damageable  Fibrous  Thermoviscoelastic  Composites. 
First  Quaterly  Report  on  SPC  96-4012.  Brussels,  1994. 

2.  N.N. Smirnov,  A. B. Kiselev,  I.D.Dimitrienko,  V.R.Dushin,  V.F. Nikitin.  Mathema¬ 
tical  modelling  of  dynamical  deforming  and  fracture  of  two-phase  composites  with 
laminated  and  fibrous  structure.  Second  Quaterly  Report  on  SPC  96-4012.  Brussels, 
1994. 

3.  N.N. Smirnov,  A. B. Kiselev,  I.D.Dimitrienko,  V.R.Dushin,  V.F. Nikitin.  Experimen¬ 
tal  method  of  determining  rheological  parameters  of  the  model  for  damaged  thermo¬ 
viscoelastic  two-phase  composites  with  laminated  structure.  Third  Quaterly  Report 
on  SPC  96-4012.  Brussels,  1994. 

4.  O.Pironeau,  B.Mohammadi.  Analysis  of  the  K-Epsilon  turbulence  model.  Mason 
Editeur,  Paris,  1994. 

5.  D. A. Anderson,  J.C.Tannehill,  R.H.Pletcher.  Computational  Fluid  Mechanics  and 
Heat  Transfer.  Hemisphere  Publ.Co.,  1984. 

6.  E.S.Oran,  J.P.Boris.  Numerical  simulation  of  reactive  flow.  Elsevier  Publ.Co., New 
York- Amsterdam-  London,  1987. 

7.  A. R. Mitchell,  R.Wait.  The  finite  element  method  in  partial  differential  equations. 
John  Wiley  &  Sons  Inc.,  Chichester-New  York-Brisbane- Toronto,  1977. 


163 


CONCLUSIONS 


The  model  of  damageable  thermoviscoelastic  composite  material  was  worked  out  during 
the  contract  period  incorporating  the  thermodynamic  breakup  criterion  based  on  the 
critical  value  of  the  total  dissipation  being  the  sum  of  mechanical  dissipation,  thermal 
dissipation  and  dissipation  in  continuous  destruction  causing  accumulation  of  damages  in 
tension,  shear  and  delamination.  The  detailed  description  of  the  general  mathematical 
model  and  methods  to  determine  the  critical  dissipation  can  be  found  in  the  interim 
reports. 

The  present  report  contains  the  results  of  numerical  investigations  of  behaviour  of 
a  composite  shell  of  a  containment  under  the  influence  of  internal  dynamical  loading 
incorporating  the  worked  out  model  for  a  damageable  laminate.  To  demonstrate  the 
utility  of  the  mathematical  model  it  was  incorporated  into  a  hydrocode  and  applied  for 
the  solution  of  the  problem  of  dynamical  deforming  of  a  composite  shell  in  a  uniform  and 
nonuniform  loading  caused  by  the  internal  explosion. 

The  model  makes  it  possible  to  determine  the  formation  and  growth  of  potentially 
damaged  zones  taking  into  account  nonuniformities  of  shell’s  parameters  distribution  due 
to  loading  and/or  curvature  nonuniformities.  The  material  constants  for  the  model  had 
been  taken  arbitrary.  Thus  the  results  of  the  described  numerical  experiments  are  mostly 
qualitative.  Analysis  of  the  results  shows  the  utility  of  the  worked  out  mathematical 
model. 
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Problem  parameters  file 

Tcihlc 

Thick  membrane  of  composite  materials  under  loading 

50 

Gas:  grid  points  along  OX 

50 

Gas;  grid  points  along  OR 

3.00000e+000 

Gas:  total  domain  length 

l.OOOOOe+000 

Gas:  total  domain  radius 

l.OOOOOe+005 

Gas:  initial  pressure 

3.00000e+002 

Gas:  initial  temperature 

2.90000e-002 

Gas;  molar  weight 

1.40000e+000 

Gas;  adhiabatic  ratio 

l.OOOOOe+007 

Pressure  behavior  constants:  amplitude 

l.OOOOOe+001 

charact.radius 

l.OOOOOe-004 

charact.time 

601 

Wall:  number  of  nodes 

l.OOOOOe+000 

Radius 

2.00000e+000 

Main  section  length 

2.00000e-001 

Corners'  radius 

l.OOOOOe-002 

Membrane  thickness 

6,10000e-001 

Volume  fraction  of  component  1 

9.40000e+010 

Compl;  Lame  coefficient  La 

5.30000e+010 

Compl:  Lame  coefficient  Mu 

2.00000e-005 

Compl:  Volume  thermal  extendibility 

l.OOOOOe+003 

Compl:  Spec,  heat  capacity  at  constant  deformations 

1.60000e+003 

Compl;  Density 

l.OOOOOe-004 

Compl:  Relaxation  time 

7.54000e+009 

Comp2;  Lame  coefficient  La 

1.89000e+009 

Comp2;  Lame  coefficient  Mu 

1.50000e-004 

Comp2;  Volume  thermal  extendibility 

2.00000e+003 

Comp2:  Spec,  heat  capacity  at  constant  deformations 

1.30000e+003 

Comp2:  Density 

l.OOOOOe-002 

Comp2:  Relaxation  time 

2.73000e+002 

Initial  temperature 

l.OOOOOe+000 

Unknown  damadge  constant  Big  Lambda 

l.OOOOOe+000 

Unknown  damadge  constant  Big  Lambda-Delta 

l.OOOOOe+004 

Unknown  damadge  constant  Big  Omega 

l.OOOOOe+000 

Unknown  damadge  constant  Big  A 

l.OOOOOe+004 

Unknown  damadge  constant  Big  C 

l.OOOOOe+004 

Unknown  damadge  constant  Big  D 

1.00000e-002 

Compl:  damadge  constant  E_* 

7.00000e-003 

Compl:  damadge  constant  E^t_* 

3.00000e-003 

Comp2;  damadge  constant  E_* 

2.00000e-003 

Comp2:  damadge  constant  E'"'t_* 

5.00000e-003 

Damadge  constant  Delta_* 

l.OOOOOe+007 

Maximal  dissipated  energy 

11296 
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